Unidad Práctica 4: Visualización de Mapas¶
En este notebook trabajaremos con información geográfica. Seguiremos estudiando el data set de viajes de la encuesta origen-destino, esta vez enfocándonos en distintos patrones geográficos que puedan ayudarnos a responder preguntas específicas.
Preámbulo y Carga de Datos¶
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
import geopandas as gpd
# esto configura la calidad de la imagen. dependerá de tu resolución. el valor por omisión es 80
mpl.rcParams["figure.dpi"] = 96
# esto depende de las fuentes que tengas instaladas en el sistema.
mpl.rcParams["font.family"] = "Fira Sans Extra Condensed"
# Importar la biblioteca seaborn para la visualización de datos
import seaborn as sns
# Establecer el estilo de fondo de las gráficas como "whitegrid" en seaborn.
sns.set_style("whitegrid")
def decode_column(
df,
fname,
col_name,
index_col="Id",
value_col=None,
sep=";",
encoding="utf-8",
index_dtype=np.float64,
):
"""
Agrega una columna extra al DataFrame `df` decodificando valores desde un archivo externo.
:param df: DataFrame al que se agregará la columna extra.
:param fname: Nombre del archivo que contiene los valores a decodificar.
:param col_name: Nombre de la columna en el DataFrame `df` que queremos decodificar.
:param index_col: Nombre de la columna en el archivo `fname` que contiene los índices que codifican la columna `col_name`.
:param value_col: Nombre de la columna en el archivo `fname` que tiene los valores decodificados. Si no se proporciona, se utilizará "value" por defecto.
:param sep: Carácter que separa los valores en el archivo `fname`.
:param encoding: Identificación del conjunto de caracteres (character set) que utiliza el archivo. Usualmente es utf-8, si no funciona, se puede probar con iso-8859-1.
:param index_dtype: Tipo de datos del índice en el archivo `fname`, por defecto np.float64.
:return: Serie con los valores decodificados que se agregarán como columna extra en el DataFrame `df`.
"""
if value_col is None:
value_col = "value"
# Lee el archivo `fname` que contiene los valores decodificados junto con sus índices.
values_df = pd.read_csv(
fname,
sep=sep,
index_col=index_col,
names=[index_col, value_col],
header=0,
dtype={index_col: index_dtype},
encoding=encoding,
)
# Extrae la columna `col_name` del DataFrame `df`.
src_df = df.loc[:, (col_name,)]
# Une el DataFrame original con los valores decodificados usando la columna `col_name` como índice.
# Se obtiene una Serie con los valores decodificados que se agregarán como columna extra en el DataFrame `df`.
return src_df.join(values_df, on=col_name)[value_col]
# leer personas
path = "https://raw.githubusercontent.com/zorzalerrante/aves/master/data/external/EOD_STGO/"
personas = pd.read_csv(
path + "personas.csv",sep=";", decimal=",", encoding="utf-8"
).rename(columns={"Factor": "FactorPersona"})
# agregar atributos
atributos = ["Sexo","TramoIngreso","Relacion","Ocupacion"]
for col_name in atributos:
personas[col_name] = decode_column(personas,path + f"Tablas_parametros/{col_name}.csv", col_name)
# leer viajes
viajes = (
pd.read_csv(path + "viajes.csv", sep=";", decimal=",")
.join(
pd.read_csv(path + "ViajesDifusion.csv", sep=";", index_col="Viaje"),
on="Viaje",
)
.join(
pd.read_csv(path + "DistanciaViaje.csv", sep=";", index_col="Viaje"),
on="Viaje",
)
)
# agregar atributos
viajes["ModoAgregado"] = decode_column(
viajes,path + "Tablas_parametros/ModoAgregado.csv", index_col="ID",value_col="Modo",col_name = "ModoAgregado")
viajes["ModoDifusion"] = decode_column(
viajes,path + "Tablas_parametros/ModoDifusion.csv", encoding="latin-1",index_col="ID",col_name = "ModoDifusion")
viajes["SectorOrigen"] = decode_column(
viajes,path + "Tablas_parametros/Sector.csv", col_name="SectorOrigen",index_col="Sector",value_col="Nombre",sep=";")
viajes["SectorDestino"] = decode_column(
viajes,path + "Tablas_parametros/Sector.csv", col_name="SectorDestino",index_col="Sector",value_col="Nombre",sep=";")
viajes["Proposito"] = decode_column(
viajes,path + "Tablas_parametros/Proposito.csv", col_name="Proposito")
viajes["ComunaOrigen"] = decode_column(
viajes,path + "Tablas_parametros/Comunas.csv","ComunaOrigen",value_col="Comuna",sep=",")
viajes["ComunaDestino"] = decode_column(
viajes,path + "Tablas_parametros/Comunas.csv", "ComunaDestino",value_col="Comuna",sep=",")
viajes["Periodo"] = decode_column(
viajes,path + "Tablas_parametros/Periodo.csv","Periodo",sep=";",value_col="Periodos")
# correciones
viajes = viajes[pd.notnull(viajes["HoraIni"])]
viajes = viajes[viajes["Imputada"] == 0].copy()
viajes["HoraIni"] = pd.to_timedelta(viajes["HoraIni"] + ":00")
# Juntar Datos
viajes_persona = viajes.merge(personas)
viajes_persona["PesoLaboral"] = (
viajes_persona["FactorLaboralNormal"] * viajes_persona["Factor_LaboralNormal"]
)
viajes_persona = viajes_persona[pd.notnull(viajes_persona["PesoLaboral"])]
print(
"{} viajes expandidos a {}".format(
len(viajes_persona), int(viajes_persona["PesoLaboral"].sum())
)
)
65591 viajes expandidos a 14566372
¿Dónde se concentran las personas que utilizan cada modo de transporte en la ciudad?¶
Es de interés saber dónde viven las personas que utilizan cada modo de transporte. Eso permite informar la planificación de nuevas redes de transporte y la gestión de las redes actuales.
Esta pregunta tiene una componente geográfica: el dónde. Si tuviésemos un enfoque de tablas, podríamos utilizar el atributo Comuna o Sector y resolver tareas del tipo parte-de-un-todo para conocer la distribución del uso de modo de transporte a través de las comunas o sectores.
Sin embargo, existe variabilidad dentro de les habitantes de cada comuna. Además, también es de interés saber si dos comunas/sectores/barrios/etc. cercanos tienen comportamiento similar, y si no, entender el contexto urbano que causa la diferencia.
Ante esa necesidad, un mapa es inmejorable. Pero, ¿cómo construir y configurar el mapa? Exploremos eso ahora.
Primero, veamos cómo luces las coordenadas que vienen dentro de los datos:
viajes_persona[
["OrigenCoordX", "OrigenCoordY", "DestinoCoordX", "DestinoCoordY"]
].head()
| OrigenCoordX | OrigenCoordY | DestinoCoordX | DestinoCoordY | |
|---|---|---|---|---|
| 0 | 335208.7188 | 6288387.0 | 338812.3125 | 6292391.0 |
| 1 | 338812.2813 | 6292391.0 | 335208.7188 | 6288387.0 |
| 2 | 338536.4375 | 6291928.0 | 354267.3438 | 6302297.0 |
| 3 | 354267.3438 | 6302297.0 | 338536.4375 | 6291928.0 |
| 4 | 338536.4375 | 6291928.0 | 350841.6563 | 6297212.0 |
Si bien las coordenadas están en formato numérico, no las tenemos estructuradas en un GeoDataFrame. Para ello, usaremos la función to_point_geodataframe para darle contexto geográfico a los datos. Crearemos dos estructuras, una para los orígenes de los viajes y otra para los destinos. Del notebook anterior ya conocemos el sistema de coordenadas:
import shapely
def to_point_geodataframe(df, longitude, latitude, crs="epsg:4326", drop=False):
gdf = gpd.GeoDataFrame(
df,
geometry=df.apply(
lambda x: shapely.geometry.Point((x[longitude], x[latitude])), axis=1
),
crs=crs,
)
if drop:
gdf = gdf.drop([longitude, latitude], axis=1)
return gdf
origenes_viajes = to_point_geodataframe(
viajes_persona, "OrigenCoordX", "OrigenCoordY", crs="epsg:32719"
)
destinos_viajes = to_point_geodataframe(
viajes_persona, "DestinoCoordX", "DestinoCoordY", crs="epsg:32719"
)
Cargamos el archivo GeoJSON que creamos en el notebook 03-python-mapas-preliminario.ipynb:
path = "https://raw.githubusercontent.com/zorzalerrante/aves/master/data/processed/scl_zonas_urbanas.json"
zones = (
gpd.read_file(path)
.set_index("ID")
.to_crs(origenes_viajes.crs)
)
zones.head()
| AREA | Zona | Com | Comuna | REGION | NOM_REGION | PROVINCIA | NOM_PROVIN | NOM_COMUNA | URBANO | TIPO | NOM_CATEG | SHAPE_Leng | SHAPE_Area | area_m2 | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ID | ||||||||||||||||
| 103 | 837.7500 | 103.0 | 13105 | El Bosque | 13 | REGIÓN METROPOLITANA DE SANTIAGO | 131 | SANTIAGO | EL BOSQUE | EL BOSQUE | CAPITAL COMUNAL | CIUDAD | 0.152123 | 0.001391 | 4.575649e+05 | POLYGON ((346371.218 6285956.953, 346349.534 6... |
| 104 | 998.8125 | 104.0 | 13105 | El Bosque | 13 | REGIÓN METROPOLITANA DE SANTIAGO | 131 | SANTIAGO | EL BOSQUE | EL BOSQUE | CAPITAL COMUNAL | CIUDAD | 0.152123 | 0.001391 | 7.725462e+05 | POLYGON ((344415.434 6285878.048, 344548.462 6... |
| 106 | 286.2500 | 106.0 | 13105 | El Bosque | 13 | REGIÓN METROPOLITANA DE SANTIAGO | 131 | SANTIAGO | EL BOSQUE | EL BOSQUE | CAPITAL COMUNAL | CIUDAD | 0.152123 | 0.001391 | 2.693838e+06 | POLYGON ((342466.924 6284938.448, 342467.002 6... |
| 115 | 857.4375 | 115.0 | 13105 | El Bosque | 13 | REGIÓN METROPOLITANA DE SANTIAGO | 131 | SANTIAGO | EL BOSQUE | EL BOSQUE | CAPITAL COMUNAL | CIUDAD | 0.152123 | 0.001391 | 7.534193e+05 | POLYGON ((344391.421 6285910.048, 344365.314 6... |
| 116 | 853.9375 | 116.0 | 13105 | El Bosque | 13 | REGIÓN METROPOLITANA DE SANTIAGO | 131 | SANTIAGO | EL BOSQUE | EL BOSQUE | CAPITAL COMUNAL | CIUDAD | 0.152123 | 0.001391 | 7.184305e+05 | POLYGON ((345371.399 6286207.024, 345608.811 6... |
Ahora nos aseguramos de tener orígenes y destinos que solo se dan dentro del contexto urbano determinado por el GeoDataFrame de zonas:
def clip_point_geodataframe(geodf, bounding_box, buffer=0):
bounds = shapely.geometry.box(*bounding_box).buffer(buffer)
return geodf[geodf.within(bounds)]
origenes_viajes = origenes_viajes[
origenes_viajes["Viaje"].isin(destinos_viajes["Viaje"])
]
origenes_viajes = clip_point_geodataframe(origenes_viajes, zones.total_bounds)
destinos_viajes = destinos_viajes[
destinos_viajes["Viaje"].isin(origenes_viajes["Viaje"])
]
Grafiquemos los datos. Utilizaremos la función figure_from_geodataframe de aves que configura los ejes de manera automática, y el método plot de geopandas puede graficar directamente los puntos de origen/destino:
def figure_from_geodataframe(
geodf,
height=5,
bbox=None,
remove_axes=True,
set_limits=True,
basemap=None,
basemap_interpolation="hanning",
):
if bbox is None:
bbox = geodf.total_bounds
aspect = (bbox[2] - bbox[0]) / (bbox[3] - bbox[1])
fig, ax = plt.subplots(figsize=(height * aspect, height))
if set_limits:
ax.set_xlim([bbox[0], bbox[2]])
ax.set_ylim([bbox[1], bbox[3]])
# code from geopandas
if geodf.crs and geodf.crs.is_geographic:
bounds = geodf.total_bounds
y_coord = np.mean([bounds[1], bounds[3]])
ax.set_aspect(1 / np.cos(y_coord * np.pi / 180))
# formula ported from R package sp
# https://github.com/edzer/sp/blob/master/R/mapasp.R
else:
ax.set_aspect("equal")
if remove_axes:
ax.set_axis_off()
if basemap is not None:
cx.add_basemap(
ax,
crs=geodf.crs.to_string(),
source=basemap,
interpolation=basemap_interpolation,
zorder=0,
)
return fig, ax
fig, ax = figure_from_geodataframe(zones, height=9, remove_axes=True)
zones.plot(ax=ax, color="#efefef", edgecolor="#abacab", linewidth=1)
origenes_viajes.plot(ax=ax, markersize=1, marker=".", color="blue", alpha=0.05)
destinos_viajes.plot(ax=ax, markersize=1, marker=".", color="red", alpha=0.05)
plt.show()
Nos damos cuenta que hay muchos datos. Tenemos que:
- Filtrarlos, es decir, decidir cuáles son relevantes para la tarea.
- Visualizarlos utilizando la técnica adecuada.
Para filtrarlos, fíjemonos en cuatro modos de transporte: público (Bip!), privado (Auto) y activo (Caminata y Bicicleta). Y por ahora veamos los viajes al trabajo.
Para visualizarlos, en aves tenemos una función dot_map que hace un mapa de puntos y que nos permitirá configurar algunos aspectos del gráfico considerando los aspectos que vimos en la clase teórica.
El gráfico se puede hacer así:
def bubble_map(
ax,
geodf: gpd.GeoDataFrame,
size,
scale=1,
palette=None,
color=None,
add_legend=True,
edgecolor="white",
alpha=1.0,
label=None,
**kwargs
):
marker = "o"
if size is not None:
if type(size) == str:
marker_size = geodf[size]
else:
marker_size = float(size)
else:
marker_size = 1
return geodf.plot(
ax=ax,
marker=marker,
markersize=marker_size * scale,
edgecolor=edgecolor,
alpha=alpha,
facecolor=color,
legend=add_legend,
label=label,
**kwargs
)
def dot_map(
ax,
geodf: gpd.GeoDataFrame,
size=10,
palette=None,
add_legend=True,
label=None,
**kwargs
):
return bubble_map(
ax,
geodf,
size=float(size),
palette=palette,
add_legend=add_legend,
edgecolor="none",
label=label,
**kwargs
)
fig, ax = figure_from_geodataframe(zones, height=6, remove_axes=True)
zones.plot(ax=ax, color="#efefef", edgecolor="#abacab", linewidth=1)
for idx, group in origenes_viajes[(origenes_viajes.Proposito == "Al trabajo")].groupby(
"ModoDifusion"
):
if idx in ["Bip!", "Auto", "Caminata", "Bicicleta"]:
dot_map(ax, group, size=10, label=idx)
ax.legend()
plt.show()
Notamos que cada categoría utiliza un tono de color distinto. Sin embargo, cuesta ver una posible distribución geográfica debido a la oclusión y aglomeramiento de los viajes.
Una posible solución es darle un tamaño a cada punto que sea proporcional al peso de los viajes. En tal caso, en vez de dot_map utilizaremos bubble_map:
fig, ax = figure_from_geodataframe(zones, height=9, remove_axes=True)
zones.plot(ax=ax, color="#efefef", edgecolor="white", linewidth=1)
for idx, group in origenes_viajes[(origenes_viajes.Proposito == "Al trabajo")].groupby(
"ModoDifusion"
):
if idx in ["Bip!", "Auto", "Caminata", "Bicicleta"]:
bubble_map(ax, group, size="PesoLaboral", scale=0.15, label=idx)
ax.legend()
fig.tight_layout()
plt.show()
La idea sonaba bien, sin embargo, el gráfico no nos permite responder la pregunta (aunque es un lindo gráfico para compartir y contemplar). Esto se debe a varios factores, incluyendo la interferencia entre los canales utilizados en la codificación.
También la implementación del gráfico es compleja. Requiere ciclos a través de los resultados de una operación groupby.
Utilizaremos la grilla GeoFacetGrid de aves (inspirada en FacetGrid de seaborn) para hacer el mismo gráfico pero utilizando conceptos vistos en clase.
from matplotlib_scalebar.scalebar import (
ScaleBar,
)
def north_arrow(
ax,
x=0.98,
y=0.06,
arrow_length=0.04,
text="N",
font_name=None,
font_size=None,
color="#000000",
arrow_color="#000000",
arrow_width=3,
arrow_headwidth=7,
):
ax.annotate(
text,
xy=(x, y),
xytext=(x, y - arrow_length),
arrowprops=dict(
facecolor=arrow_color, width=arrow_width, headwidth=arrow_headwidth
),
ha="center",
va="center",
fontsize=font_size,
fontname=font_name,
color=color,
xycoords=ax.transAxes,
)
def geographical_scale(ax, location="lower left"):
ax.add_artist(ScaleBar(1, location="lower left"))
from seaborn.axisgrid import FacetGrid, Grid
from matplotlib import colorbar
class GeoFacetGrid(FacetGrid):
def __init__(self, geodataframe: gpd.GeoDataFrame, *args, **kwargs):
geocontext = kwargs.pop("context", None)
if geocontext is None:
geocontext = geodataframe
self.geocontext = geocontext
self.bounds = geocontext.total_bounds
self.aspect = (self.bounds[2] - self.bounds[0]) / (
self.bounds[3] - self.bounds[1]
)
kwargs["aspect"] = self.aspect
kwargs["xlim"] = (self.bounds[0], self.bounds[2])
kwargs["ylim"] = (self.bounds[1], self.bounds[3])
super().__init__(geodataframe, *args, **kwargs)
for ax in self.axes.flatten():
if kwargs.get("remove_axes", True):
ax.set_axis_off()
aspect = kwargs.get("aspect", "auto")
# code from geopandas
if aspect == "auto":
if geodataframe.crs and geodataframe.crs.is_geographic:
bounds = geodataframe.total_bounds
y_coord = np.mean([bounds[1], bounds[3]])
ax.set_aspect(1 / np.cos(y_coord * np.pi / 180))
# formula ported from R package sp
# https://github.com/edzer/sp/blob/master/R/mapasp.R
else:
ax.set_aspect("equal")
elif aspect is not None:
ax.set_aspect(aspect)
self.zorder = 0
def add_layer(self, func_or_data, *args, **kwargs):
if isinstance(func_or_data, gpd.GeoDataFrame):
# a direct geography
for ax in self.axes.flatten():
func_or_data.plot(*args, ax=ax, zorder=self.zorder, **kwargs)
else:
plot = lambda *a, **kw: func_or_data(
plt.gca(), kw.pop("data"), *a, zorder=self.zorder, **kw
)
self.map_dataframe(plot, *args, **kwargs)
self.zorder += 1
def add_basemap(
self, file_path, interpolation="hanning", reset_extent=False, **kwargs
):
for ax in self.axes.flatten():
cx.add_basemap(
ax,
crs=self.geocontext.crs.to_string(),
source=file_path,
interpolation=interpolation,
zorder=self.zorder,
reset_extent=reset_extent,
**kwargs
)
if not reset_extent:
ax.set_xlim((self.bounds[0], self.bounds[2]))
ax.set_ylim((self.bounds[1], self.bounds[3]))
self.zorder += 1
def add_map_elements(
self,
all_axes=False,
scale=True,
arrow=True,
scale_args={},
arrow_args={},
):
for ax in self.axes.flatten():
if arrow:
north_arrow(ax, **arrow_args)
if scale:
geographical_scale(ax, **scale_args)
if not all_axes:
break
def add_global_colorbar(self, palette, k, title=None, title_args={}, **kwargs):
orientation = kwargs.get("orientation", "horizontal")
if orientation == "horizontal":
cax = self.fig.add_axes([0.25, -0.012, 0.5, 0.01])
elif orientation == "vertical":
cax = self.fig.add_axes([1.01, 0.25, 0.01, 0.5])
else:
raise ValueError("unsupported colorbar orientation")
if title:
cax.set_title(title, **title_args)
cb = colorbar.ColorbarBase(
cax, cmap=colormap_from_palette(palette, n_colors=k), **kwargs
)
cax.set_axis_off()
return cax, cb
def set_title(self, title, **kwargs):
self.fig.suptitle(title, **kwargs)
#from aves.visualization.figures import GeoFacetGrid
grid = GeoFacetGrid(
# los datos
origenes_viajes,
# el contexto geográfico. es opcional, se utiliza para configurar el gráfico
context=zones,
# la variable que mapearemos a las filas del gráfico
row="Proposito",
# en este caso, solo una fila
row_order=["Al trabajo"],
# el canal de codificación tono (hue) expresará la columna categórica ModoDifusion
hue="ModoDifusion",
# solo veremos estas categorías
hue_order=["Auto", "Bicicleta", "Bip!", "Caminata"],
# la altura del gráfico
height=9,
)
# visualizamos el contexto
grid.add_layer(zones, color="#efefef", edgecolor="white", linewidth=1)
# agregamos los bubble_map correspondientes
grid.add_layer(bubble_map, size="PesoLaboral", scale=0.15)
# agregamos la leyenda
grid.add_legend()
grid.fig.tight_layout()
La implementación del gráfico es distinta, sin embargo, los problemas conceptuales se mantienen. Ahora bien, al utilizar la grilla, podemos utilizar sus opciones para reducir su complejidad a través de facetas.
Por ejemplo, podemos dividir el gráfico en columnas, donde cada columna muestra un modo de transporte diferente:
grid = GeoFacetGrid(
origenes_viajes,
context=zones,
row="Proposito",
col="ModoDifusion",
row_order=["Al trabajo", "De salud"],
col_order=["Auto", "Bip!", "Caminata", "Bicicleta"],
height=9,
hue="ModoDifusion",
palette="plasma",
)
grid.add_layer(zones, color="#efefef", edgecolor="white", linewidth=1)
grid.add_layer(bubble_map, size="PesoLaboral", scale=0.15, edgecolor="black", alpha=0.5)
¡Mucho mejor! Sin embargo, la sobreposición de las burbujas sigue siendo un problema. No logramos ver la distribución debido a la oclusión de los puntos.
La solución es usar un heat_map. Así se ve un heatmap sobre todos los viajes al trabajo. Solo debemos reemplazar el uso de bubble_map por heat_map (y usar los parámetros correspondientes):
import KDEpy
import matplotlib.colors as colors
def kde_from_points(
geodf,
kernel="gaussian",
norm=2,
bandwidth=1e-2,
grid_points=2 ** 9,
weight_column=None,
):
# La variable grid_points define la cantidad de puntos en el espacio en el que se estimará la densidad
# hacemos una lista con las coordenadas de los viajes
point_coords = np.vstack([geodf.geometry.x, geodf.geometry.y]).T
# instanciamos la Fast-Fourier Transform Kernel Density Estimation
kde = KDEpy.FFTKDE(bw=bandwidth, norm=norm, kernel=kernel)
weights = None if weight_column is None else geodf[weight_column].values
grid, points = kde.fit(point_coords, weights=weights).evaluate(grid_points)
x, y = np.unique(grid[:, 0]), np.unique(grid[:, 1])
z = points.reshape(grid_points, grid_points).T
return x, y, z
def colormap_from_palette(palette_name, n_colors=10):
return colors.ListedColormap(sns.color_palette(palette_name, n_colors=n_colors))
def heat_map(
ax,
geodf,
weight=None,
low_threshold=0,
max_threshold=1.0,
n_levels=5,
alpha=1.0,
palette="magma",
kernel="cosine",
norm=2,
bandwidth=1e-2,
grid_points=2 ** 6,
return_heat=False,
cbar_label=None,
cbar_width=2.4,
cbar_height=0.15,
cbar_location="upper left",
cbar_orientation="horizontal",
cbar_pad=0.05,
cbar_bbox_to_anchor=(0.0, 0.0, 1.0, 1.0),
cbar_bbox_transform=None,
legend_type="none",
**kwargs
):
heat = kde_from_points(
geodf,
kernel=kernel,
norm=norm,
bandwidth=bandwidth,
grid_points=grid_points,
weight_column=weight,
)
norm_heat = heat[2] / heat[2].max()
cmap = colormap_from_palette(palette, n_colors=n_levels)
levels = np.linspace(low_threshold, max_threshold, n_levels)
# TODO: this should be factorized into an utility function
if legend_type == "colorbar":
# add_ranged_color_legend(ax)
if cbar_bbox_transform is None:
cbar_bbox_transform = ax.transAxes
if cbar_location != "out":
cbar_ax = inset_axes(
ax,
width=cbar_width,
height=cbar_height,
loc=cbar_location,
bbox_to_anchor=cbar_bbox_to_anchor,
bbox_transform=cbar_bbox_transform,
borderpad=0,
)
else:
divider = make_axes_locatable(ax)
cbar_main = divider.append_axes(
"bottom" if cbar_orientation == "horizontal" else "right",
size=cbar_height if cbar_orientation == "horizontal" else cbar_width,
pad=cbar_pad,
)
cbar_main.set_axis_off()
cbar_ax = inset_axes(
cbar_main,
width=cbar_width,
height=cbar_height,
loc="center",
bbox_to_anchor=(0.0, 0.0, 1.0, 1.0),
bbox_transform=cbar_main.transAxes,
borderpad=0,
)
color_legend(cbar_ax, cmap, levels, orientation=cbar_orientation)
else:
cbar_ax = None
if not return_heat:
return (
ax.contourf(heat[0], heat[1], norm_heat, levels, alpha=alpha, cmap=cmap),
cbar_ax,
)
else:
return (
ax.contourf(heat[0], heat[1], norm_heat, levels, alpha=alpha, cmap=cmap),
cbar_ax,
heat,
)
grid = GeoFacetGrid(
origenes_viajes,
context=zones,
row="Proposito",
col="ModoDifusion",
row_order=["Al trabajo"],
col_order=["Auto", "Bip!", "Caminata", "Bicicleta"],
height=9,
hue="ModoDifusion"
)
grid.add_layer(zones, color="#efefef", edgecolor="white", linewidth=1)
grid.add_layer(
heat_map,
# atributo de los datos con la importancia o peso de cada viaje
weight="PesoLaboral",
# cantidad de niveles/colores del mapa de calor
n_levels=10,
# radio de influencia de cada viaje
bandwidth=1000,
# valor de corte para los valores bajos del heatmap
low_threshold=0.025,
# transparencia
alpha=0.75,
# paleta de colores
palette="inferno"
)
Esta técnica de visualización responde nuestra pregunta: nos dice dónde hay concentración, y, al mismo tiempo, deduce una forma, unos límites, para esos lugares que concentran actividad.
Pero hay algo extraño en el gráfico. Solo vemos concentración en uno o dos lugares de la ciudad. Queríamos analizar de acuerdo al lugar de residencia de las personas, en un tipo de viaje importante, pero no estamos viendo mapas que nos hagan pensar en residencia. ¿Qué sucede?
La respuesta es que nos fijamos en los orígenes de los viajes al trabajo, pero sin discriminar que el viaje al trabajo fuese desde el hogar de las personas. Y, en un día laboral normal, usualmente se hacen dos viajes al trabajo. Uno desde la casa, y otro al regresar al trabajo luego de comer.
Hay soluciones para esto. La primera sería corregir ese análisis. No lo haremos puesto que no es una clase de análisis de datos. Queda como problema propuesto.
La segunda es utilizar los destinos de los viajes. En un día laboral normal podemos asumir que es común que solo se regrese a casa una vez, y que ese regreso se realice en el modo de transporte principal.
grid = GeoFacetGrid(
destinos_viajes,
context=zones,
col="ModoDifusion",
row="Proposito",
row_order=["volver a casa"],
col_order=["Auto", "Bip!", "Caminata", "Bicicleta"],
height=9,
)
grid.add_layer(zones, color="#efefef", edgecolor="white", linewidth=1)
grid.add_layer(
heat_map,
weight="PesoLaboral",
n_levels=10,
bandwidth=1000,
low_threshold=0.025,
alpha=0.75,
palette="inferno",
)
Configuremos el gráfico para que sea más amigable. Haremos lo siguiente:
- Haremos un gráfico de 2x2. Para ello necesitamos filtrar por nuestra cuenta los datos, de modo que el objeto
gridpueda disponer de las filas (o de las columnas) repitiendo valores. - Agregaremos como base la imagen del territorio que descargamos en el notebook anterior.
- Utilizaremos una leyenda de colores para explicar cuál es el color asociado a la concentración.
import contextily as cx
grid = GeoFacetGrid(
destinos_viajes[destinos_viajes["Proposito"] == "volver a casa"],
context=zones,
col="ModoDifusion",
col_wrap=2,
col_order=["Auto", "Bip!", "Caminata", "Bicicleta"],
height=9,
)
grid.add_basemap("data/scl_toner_12.tif")
grid.add_layer(
heat_map,
weight="PesoLaboral",
n_levels=10,
bandwidth=1000,
low_threshold=0.025,
alpha=0.75,
palette="inferno",
)
grid.add_global_colorbar(
"inferno",
10,
title="Intensidad de Viajes (de menos a más)",
orientation="horizontal",
)
grid.tight_layout()
# grid.savefig('../reports/figures/example_geofacetgrid.png', dpi=150)
C:\Users\franc\AppData\Local\pypoetry\Cache\virtualenvs\db-connectors-RsEieBu8-py3.8\lib\site-packages\seaborn\axisgrid.py:118: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect. self._figure.tight_layout(*args, **kwargs)
<__main__.GeoFacetGrid at 0x18c0e988550>
Agregar el contexto territorial nos permite ver cosas como el rol segregador que juega la autopista Américo Vespucio. ¿Qué otras cosas ven ustedes?
Ejercicio propuesto:
- Consideren el análisis por grupos etáreos que hicimos la clase pasada. Definan categorías de grupos y determinen cuáles lugares de la ciudad visitan, tanto en general como por propósito de viaje. Rangos de ejemplo incluyen: menores de edad (< 18 años), adultos mayores (> 65 años), etc.
Nota: Debes tener descargado el archivo scl_toner_12.tif en tu computador local o google colab.
¿Cuán lejos queda el trabajo de acuerdo al lugar de residencia?¶
Con esta pregunta queremos entender si existe un patrón geográfico en las elecciones de residencia y trabajo de las personas.
Para responder la pregunta, primero filtramos los viajes que nos interesan:
viajes_trabajo = origenes_viajes[(origenes_viajes.Proposito == 'Al trabajo') &
(pd.notnull(origenes_viajes.PesoLaboral)) &
(origenes_viajes.DistEuclidiana > 0)].drop_duplicates(subset='Persona', keep='first')
print(len(viajes_trabajo), viajes_trabajo.PesoLaboral.sum())
11238 2380126.1986132185
Observamos que, si bien son 11228 filas en la tabla.
La columna DistEuclidiana contiene la distancia entre los puntos de origen y destino de los viajes.
No todos los viajes al trabajo se efectúan desde la casa, ni una persona hace un único viaje al trabajo durante el día. Como una manera de alivianar el problema nos quedamos con un único viaje al trabajo por persona.
Veamos la distribución de la distancia:
viajes_trabajo['DistEuclidiana'].describe()
count 1.123800e+04 mean 1.008620e+04 std 1.900365e+04 min 1.000000e+00 25% 3.628250e+03 50% 8.025000e+03 75% 1.414975e+04 max 1.258473e+06 Name: DistEuclidiana, dtype: float64
Ahora bien, el promedio que nos entregó el método pd.describe no considera la representatividad de cada viaje. Podemos utilizar la función weighted_mean para calcular el promedio ponderado:
def weighted_mean(df, value_column, weighs_column):
weighted_sum = (df[value_column] * df[weighs_column]).sum()
return weighted_sum / df[weighs_column].sum()
viajes_trabajo['DistEuclidiana'].mean(), weighted_mean(viajes_trabajo, 'DistEuclidiana', 'PesoLaboral')
(10086.19514148425, 9549.626344374019)
En este caso, el promedio ponderado no está tan lejos del promedio sin ponderar, pero nada asegura que esa diferencia sea pequeña en toda la ciudad.
Podemos calcular la distancia promedio al trabajo por zona haciendo una operación pd.groupby:
distancia_zonas = (viajes_trabajo
.groupby(['ZonaOrigen'])
.apply(lambda x: weighted_mean(x, 'DistEuclidiana', 'PesoLaboral'))
.rename('distancia_al_trabajo')
)
distancia_zonas
ZonaOrigen
1 6452.138966
2 5952.584070
3 4091.135703
4 4731.477408
5 17567.000000
...
857 22153.295597
858 18184.397206
859 17326.810773
860 16946.779972
861 12782.225811
Name: distancia_al_trabajo, Length: 738, dtype: float64
La serie contiene, para cada zona, la distancia promedio al trabajo de la gente que vive en ella (y posiblemente otras más). Si queremos ver esta distribución podemos utilizar el método plot(kind='kde'), como hicimos en notebooks anteriores:
distancia_zonas.plot(kind='kde')
plt.xlim([0, distancia_zonas.max()])
plt.title('Distancia al Trabajo por Zonas')
plt.xlabel('Distancia')
plt.ylabel('Densidad (KDE)')
sns.despine()
Ahora bien, esta distribución solamente nos permite saber propiedades estádisticas. ¡No tenemos un contexto geográfico! Por eso necesitamos el mapa, para saber si la distancia tiene relación con la ubicación de cada zona en la ciudad, con saber si zonas que están cerca entre sí tienen distancias al trabajo similares.
Como tenemos los valores de distancia por zona, un choropleth map nos permitiría mostrar la distancia asociada a cada zona al mismo tiempo que su posición geográfica.
zones.crs
<Projected CRS: EPSG:32719> Name: WGS 84 / UTM zone 19S Axis Info [cartesian]: - E[east]: Easting (metre) - N[north]: Northing (metre) Area of Use: - name: Between 72°W and 66°W, southern hemisphere between 80°S and equator, onshore and offshore. Argentina. Bolivia. Brazil. Chile. Colombia. Peru. - bounds: (-72.0, -80.0, -66.0, 0.0) Coordinate Operation: - name: UTM zone 19S - method: Transverse Mercator Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich
def color_legend(
ax,
color_list,
bins=None,
sizes=None,
orientation="horizontal",
remove_axes=False,
bin_spacing="proportional",
tick_labels=None,
):
if bins is None:
if type(color_list) == colors.ListedColormap:
N = color_list.N
else:
N = len(color_list)
bins = np.array(range(N))
if sizes is not None:
bar_width = bins[1:] - bins[0:-1]
if orientation == "horizontal":
ax.bar(
bins[:-1],
sizes,
width=bar_width,
align="edge",
color=color_list,
edgecolor=color_list,
)
ax.set_xticks(bins, labels=tick_labels)
else:
ax.barh(
bins[:-1],
sizes,
height=bar_width,
align="edge",
color=color_list,
edgecolor=color_list,
)
ax.set_yticks(bins, labels=tick_labels)
else:
cbar_norm = colors.BoundaryNorm(bins, len(bins) - 1)
if type(color_list) == colors.ListedColormap:
cmap = color_list
else:
cmap = colors.ListedColormap(color_list)
cb = colorbar.ColorbarBase(
ax,
cmap=cmap,
norm=cbar_norm,
ticks=bins,
spacing=bin_spacing,
orientation=orientation,
)
if tick_labels:
cb.set_ticklabels(tick_labels)
sns.despine(ax=ax, top=True, bottom=True, left=True, right=True)
if remove_axes:
ax.set_axis_off()
return ax
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
def add_ranged_color_legend(
ax,
bins,
built_palette,
location="lower center",
orientation="horizontal",
label=None,
label_size="x-small",
title_size="small",
title_align="left",
width="50%",
height="5%",
bbox_to_anchor=(0.0, 0.0, 0.95, 0.95),
bbox_transform=None,
tick_labels=None,
**kwargs,
):
if bbox_transform is None:
bbox_transform = ax.transAxes
if location != "out":
cbar_ax = inset_axes(
ax,
width=width,
height=height,
loc=location,
bbox_to_anchor=bbox_to_anchor,
bbox_transform=bbox_transform,
)
else:
divider = make_axes_locatable(ax)
cbar_main = divider.append_axes(
"bottom" if orientation == "horizontal" else "right",
size=height if orientation == "horizontal" else width,
)
cbar_main.set_axis_off()
cbar_ax = inset_axes(
cbar_main,
width=width,
height=height,
loc="center",
bbox_to_anchor=(0.0, 0.0, 1.0, 1.0),
bbox_transform=cbar_main.transAxes,
borderpad=0,
)
color_legend(
cbar_ax,
built_palette,
bins,
orientation=orientation,
remove_axes=False,
tick_labels=tick_labels if tick_labels is not None else [],
**kwargs,
)
sns.despine(ax=cbar_ax, left=True, top=True, bottom=True, right=True)
if label:
cbar_ax.set_title(label, loc="left", fontsize=title_size)
cbar_ax.tick_params(labelsize=label_size)
return cbar_ax
from mapclassify import FisherJenks, Quantiles
def choropleth_map(
ax,
geodf,
column,
k=10,
palette=None,
default_divergent="RdBu_r",
default_negative="Blues_r",
default_positive="Reds",
palette_type="light",
legend="colorbar",
edgecolor="white",
palette_center=None,
binning="uniform",
bins=None,
alpha=1.0,
linewidth=1,
zorder=1,
cbar_args={},
**kwargs,
):
geodf = geodf[pd.notnull(geodf[column])].copy()
min_value, max_value = geodf[column].min(), geodf[column].max()
if binning in ("fisher_jenks", "quantiles"):
if binning == "fisher_jenks":
binning_method = FisherJenks(geodf[column], k=k)
else:
binning_method = Quantiles(geodf[column], k=k)
bins = np.insert(binning_method.bins, 0, geodf[column].min())
geodf = geodf.assign(__bin__=binning_method.yb)
elif binning == "uniform":
bins = np.linspace(
min_value, max_value + (max_value - min_value) * 0.001, num=k + 1
)
geodf = geodf.assign(
__bin__=lambda x: pd.cut(
x[column], bins=bins, include_lowest=True, labels=False
).astype(np.int)
)
elif binning == "custom":
if bins is None:
raise ValueError("bins are needed for custom binning")
bins = np.array(bins)
geodf = geodf.assign(
__bin__=lambda x: pd.cut(
x[column], bins=bins, include_lowest=True, labels=False
).astype(np.int)
)
else:
raise ValueError(
"only fisher_jenks, quantiles and uniform binning are supported"
)
cmap_name = None
norm = None
midpoint = 0.0
using_divergent = False
built_palette = None
if palette_center is not None:
midpoint = palette_center
if palette is None:
if min_value < 0 and max_value > 0:
# divergent
cmap_name = default_divergent
norm = MidpointNormalize(vmin=min_value, vmax=max_value, midpoint=midpoint)
using_divergent = True
elif min_value < 0:
# sequential
cmap_name = default_negative
else:
# sequential
cmap_name = default_positive
else:
if not isinstance(palette, colors.Colormap):
if colors.is_color_like(palette):
if palette_type == "light":
built_palette = sns.light_palette(palette, n_colors=k)
else:
built_palette = sns.dark_palette(palette, n_colors=k)
else:
built_palette = sns.color_palette(palette, n_colors=k)
if norm is None:
if palette_center is None:
# norm = colors.Normalize(vmin=min_value, vmax=max_value)
norm = colors.BoundaryNorm(bins, k)
else:
norm = MidpointNormalize(vmin=min_value, vmax=max_value, midpoint=midpoint)
using_divergent = True
if built_palette is None:
if not using_divergent:
built_palette = sns.color_palette(cmap_name, n_colors=k)
else:
middle_idx = np.where((bins[:-1] * bins[1:]) <= 0)[0][0]
left = middle_idx
right = k - middle_idx - 1
if left == right:
built_palette = sns.color_palette(cmap_name, n_colors=k)
else:
delta = np.abs(left - right)
expanded_k = k + delta
start_idx = max(middle_idx - left, right - middle_idx)
built_palette = sns.color_palette(cmap_name, n_colors=expanded_k)[
start_idx : start_idx + k
]
for idx, group in geodf.groupby("__bin__"):
group.plot(
ax=ax,
facecolor=built_palette[idx],
linewidth=linewidth,
edgecolor=edgecolor,
alpha=alpha,
zorder=zorder,
aspect=None,
)
if legend is not None:
cbar_ax = add_ranged_color_legend(ax, bins, built_palette, **cbar_args)
else:
cbar_ax = None
return ax, cbar_ax
fig, ax = figure_from_geodataframe(zones, height=9, remove_axes=False)
ax, cax = choropleth_map(ax, zones.join(distancia_zonas, how='inner'), 'distancia_al_trabajo', k=5)
cax.set_title('Distancia al Trabajo [m]', loc='left')
ax.grid(True)
fig.tight_layout()
C:\Users\franc\AppData\Local\Temp\ipykernel_2608\3905406884.py:41: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information. Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations ).astype(np.int) C:\Users\franc\AppData\Local\Temp\ipykernel_2608\1264298429.py:7: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect. fig.tight_layout()
La leyenda explica que valores más claros tienen menores distancias al trabajo. Al mismo tiempo, observamos la concentración de los trabajos en la ciudad: las zonas más claras tienden a aparecer en Santiago, Providencia, Las Condes, y sus alrededores. Como consecuencia, las zonas periféricas tienen una distancia al trabajo mayor: casi toda la periferia está compuesta de zonas oscuras.
¿cómo se calculan los rangos que son representados por cada color?
Existen varias alternativas. La opción por omisión es dividir el intervalo de valores posibles en rangos uniformes (opción binning='uniform'). Otra es dividir las observaciones en cuantiles (opción binning='quantiles'). Y otra es utilizar el esquema llamado Natural Breaks, que busca la mejor manera de clasificar datos en un número específico de categorías (una especie de clustering K-Means en una dimensión). Lo podemos especificar en la función de dibujo con el parámetro binning='fisher_jenks' (se llama así porque lo inventó Jenks y lo mejoró Fisher).
Usemos la función small_multiples_from_geodataframe para comparar estas tres estrategias de codificación visual en el canal de color:
def small_multiples_from_geodataframe(
geodf,
n_variables,
height=5,
col_wrap=5,
bbox=None,
sharex=True,
sharey=True,
remove_axes=True,
set_limits=True,
flatten_axes=True,
aspect="auto",
basemap=None,
basemap_interpolation="hanning",
):
if n_variables <= 1:
return figure_from_geodataframe(
geodf,
height=height,
bbox=bbox,
remove_axes=remove_axes,
set_limits=set_limits,
basemap=basemap,
basemap_interpolation=basemap_interpolation,
)
if bbox is None:
bbox = geodf.total_bounds
# code from geopandas
if aspect == "auto":
if geodf.crs and geodf.crs.is_geographic:
y_coord = np.mean([bbox[1], bbox[3]])
aspect_ratio = 1 / np.cos(y_coord * np.pi / 180)
# formula ported from R package sp
# https://github.com/edzer/sp/blob/master/R/mapasp.R
else:
aspect_ratio = 1
else:
aspect_ratio = aspect
n_columns = min(col_wrap, n_variables)
n_rows = n_variables // n_columns
if n_rows * n_columns < n_variables:
n_rows += 1
fig, axes = plt.subplots(
n_rows,
n_columns,
figsize=(n_columns * height / aspect_ratio, n_rows * height),
sharex=sharex,
sharey=sharey,
squeeze=False,
)
flattened = axes.flatten()
if set_limits:
for ax in flattened:
ax.set_xlim([bbox[0], bbox[2]])
ax.set_ylim([bbox[1], bbox[3]])
for ax in flattened:
ax.set_aspect(aspect_ratio)
if remove_axes:
for ax in flattened:
ax.set_axis_off()
else:
# deactivate only unneeded axes
for i in range(n_variables, len(axes)):
flattened[i].set_axis_off()
if basemap is not None:
for ax in flattened:
cx.add_basemap(
ax,
crs=geodf.crs.to_string(),
source=basemap,
interpolation=basemap_interpolation,
zorder=0,
)
if flatten_axes:
return fig, flattened
return fig, axes
fig, axes = small_multiples_from_geodataframe(zones, 3, height=9, remove_axes=False)
for binning, ax in zip(['uniform', 'quantiles', 'fisher_jenks'], axes):
ax.set_title(f'binning = {binning}')
_, cax = choropleth_map(ax, zones.join(distancia_zonas, how='inner'), 'distancia_al_trabajo', k=5, binning=binning, linewidth=0.1)
cax.set_title('Distancia al Trabajo [m]', loc='left')
ax.grid(True)
fig.tight_layout()
C:\Users\franc\AppData\Local\Temp\ipykernel_2608\3905406884.py:41: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information. Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations ).astype(np.int) C:\Users\franc\AppData\Local\pypoetry\Cache\virtualenvs\db-connectors-RsEieBu8-py3.8\lib\site-packages\mapclassify\classifiers.py:1956: UserWarning: Numba not installed. Using slow pure python version. warnings.warn( C:\Users\franc\AppData\Local\Temp\ipykernel_2608\3215009884.py:8: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect. fig.tight_layout()
¿Cuál esquema les parece mejor (para una definición de "mejor" elegida por ustedes)? Cada uno tiene sus ventajas y desventajas.
¿Qué otro patrón observan en este mapa? ¿Qué más harían?
¿Se imaginan este mapa mezclado con los mapas anteriores? Ése es un ejercicio propuesto.
El esquema de Fisher-Jenks es utilizado con frecuencia posiblemente porque balancea los aspectos de los otros esquemas. Lo utilizaremos para la versión final del mapa, esta vez con información contextual. Además de la información territorial, agregaremos elementos típicos de los mapas cartográficos: una escala y una flecha hacia el norte.
grid = GeoFacetGrid(zones.join(distancia_zonas, how="inner"), height=9)
grid.add_basemap("data/scl_toner_12.tif")
grid.add_layer(
choropleth_map,
"distancia_al_trabajo",
k=5,
linewidth=0.5,
edgecolor="black",
binning="fisher_jenks",
palette="RdPu",
alpha=0.85,
cbar_args=dict(
label="Distancia [m]",
height="22%",
width="2%",
orientation="vertical",
location="center left",
label_size="small",
bbox_to_anchor=(0.0, 0.0, 0.9, 1.0),
),
)
grid.add_map_elements()
grid.set_title("Distancia al Trabajo")
grid.tight_layout()
plt.show()
#grid.savefig('../reports/figures/example_choropleth.png', dpi=150)
C:\Users\franc\AppData\Local\pypoetry\Cache\virtualenvs\db-connectors-RsEieBu8-py3.8\lib\site-packages\mapclassify\classifiers.py:1956: UserWarning: Numba not installed. Using slow pure python version. warnings.warn( C:\Users\franc\AppData\Local\pypoetry\Cache\virtualenvs\db-connectors-RsEieBu8-py3.8\lib\site-packages\seaborn\axisgrid.py:118: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect. self._figure.tight_layout(*args, **kwargs)
De esta manera el mapa nos permite caracterizar las zonas, y conjuntos de zonas, de acuerdo a la distancia que deben recorrer sus habitantes para ir al trabajo.
Queda propuesto desagregar este mapa por características de las personas, como su nivel socioeconómico.